Ce document montre l’essentiel des manipulations à effectuer avec R pour répondre aux questions du TP1 d’Analyse Univariée. Il suppose une connaissance préalable de la syntaxe R, de préférence à l’aide de l’IDE RStudio et les bases du package sf tutoriel ici qui implémente la norme SF comme dans PostGIS.

La connaissance du package dplyr et de sa notation à base de pipe (%>%) sera utile pour écrire des chaînes de traitements plus concises que celle listées ici.

Chargement des données.

Nous avons récupéré le fichier des données spatialisé sous le format GeoJSON. D’autres formats sont disponibles, notamment via une API.

Pour lire le GeoJSON , nous avons besoin d’installer un package spécifique : geojsonsf, qui se charge de la conversion entre des objets geographiques au format GeoJSON et leur représentation au format géographique sf, très commode.
Pour manipuler les objets spatiaux, nous devons aussi installer la librairie sf. Il nous faut ensuite charger ces librairie. Cela se fait au moyen des fonctions install.packages("nom_du_package") et library(nom_du_package)

install.packages("geojsonsf")
library(geojsonsf)
install.packages("sf")
library(sf)

Puis nous chargeons les données à l’aide de la fonction geojson_sf().

N.B. la fonction geojson_sf requiert le chemin d’accès au fichier, fourni sous sa forme relative ou absolue. Le repertoire de travail de R est défini par la commande setwd("/chemin/vers/le/fichier"). Pour connaître le repertoire de travail courant, utiliser la fonction getwd()

N.B. le nombre de lignes et les valeurs de données qui apparaissent dans ce support peuvent varier, les données étant mises à jour régulièrement suir le site opendata.paris.fr

Le chargement du fichier peut être un peu long: il contient 203818 observations. Pour connaitre la structure de l’objet , nous utilisons la fonction str() qui nous donne un aperçu du type des colonnes du dataframe que nous manipulons. C’est un dataframe particulier puisqu’il est de la classe sf : l’un de ses colonnes est dédié à la description de sa géométrie (ici, des points 2D)

arbres <- geojson_sf("les-arbres.geojson")
names(arbres) # affiche le nom des colonnes (variables)
##  [1] "remarquable"        "circonferenceencm"  "stadedeveloppement"
##  [4] "genre"              "idbase"             "arrondissement"    
##  [7] "idemplacement"      "geo_point_2d"       "geometry"          
## [10] "adresse"            "libellefrancais"    "complementadresse" 
## [13] "domanialite"        "typeemplacement"    "hauteurenm"        
## [16] "varieteoucultivar"  "espece"

Aperçu des données

La fonction head() permet d’observer les \(n\) premières lignes d’un dataframe

head(arbres, 10) 
## Simple feature collection with 10 features and 16 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 2.211229 ymin: 48.80865 xmax: 2.420644 ymax: 48.91002
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
##    remarquable circonferenceencm stadedeveloppement             genre
## 1         <NA>                53               <NA>             Tilia
## 2            0                80                  A              Acer
## 3         <NA>                 0               <NA>              Acer
## 4            0               155                  A          Aesculus
## 5            0               190                  A          Aesculus
## 6            0                30               <NA>            Betula
## 7         <NA>                 0               <NA>             Tilia
## 8            0               150                  A             Tilia
## 9         <NA>                45               <NA> x Cupressocyparis
## 10        <NA>                70               <NA>          Aesculus
##    idbase    arrondissement idemplacement
## 1  135722    PARIS 4E ARRDT      00000005
## 2  143025   PARIS 20E ARRDT  D00000029017
## 3  146409    HAUTS-DE-SEINE  A05400012006
## 4  147256   PARIS 20E ARRDT  A01500064008
## 5  148395   PARIS 17E ARRDT  A12600028008
## 6  151564   PARIS 20E ARRDT  D00000048002
## 7  168345 SEINE-SAINT-DENIS  A0760F161033
## 8  171022   PARIS 20E ARRDT  A09200059002
## 9  171785    HAUTS-DE-SEINE  A10200001004
## 10 194202   PARIS 12E ARRDT     093350004
##                          geo_point_2d                  geometry
## 1        [48.857165868,2.35768578933] POINT (2.357686 48.85717)
## 2  [48.8587364858,2.3958940307499998] POINT (2.395894 48.85874)
## 3       [48.8086477164,2.31291877994] POINT (2.312919 48.80865)
## 4  [48.863124092099997,2.39025920157] POINT (2.390259 48.86312)
## 5  [48.8977886755,2.3180730804199998] POINT (2.318073 48.89779)
## 6  [48.8628202878,2.3929664570399997] POINT (2.392966 48.86282)
## 7       [48.9053057888,2.42064407509] POINT (2.420644 48.90531)
## 8       [48.8608889982,2.38879046137]  POINT (2.38879 48.86089)
## 9       [48.9100232498,2.21122948601] POINT (2.211229 48.91002)
## 10       [48.8327722445,2.4063340346] POINT (2.406334 48.83277)
##                                                        adresse
## 1                                     JARDINIERE RUE DU TRESOR
## 2                          CIMETIERE DU PERE LACHAISE / DIV 29
## 3  CIMETIERE DE BAGNEUX / AVENUE DES ERABLES POURPRES / DIV 12
## 4      CIMETIERE DU PERE LACHAISE / AVENUE CIRCULAIRE / DIV 64
## 5         CIMETIERE DES BATIGNOLLES / AVENUE LATERALE / DIV 28
## 6                          CIMETIERE DU PERE LACHAISE / DIV 48
## 7     CIMETIERE DE PANTIN / AVENUE DES PETITS PONTS / DIV F161
## 8    CIMETIERE DU PERE LACHAISE / AVENUE DU BOULEVARD / DIV 59
## 9     CIMETIERE DE NANTERRE / AVENUE DU POURTOUR OUEST / DIV 1
## 10                                               TALUS N°33-50
##    libellefrancais complementadresse  domanialite typeemplacement
## 1          Tilleul              <NA>       Jardin           Arbre
## 2           Erable              <NA>    CIMETIERE           Arbre
## 3           Erable              <NA>    CIMETIERE           Arbre
## 4       Marronnier              <NA>    CIMETIERE           Arbre
## 5       Marronnier              <NA>    CIMETIERE           Arbre
## 6          Bouleau              <NA>    CIMETIERE           Arbre
## 7          Tilleul              <NA>    CIMETIERE           Arbre
## 8          Tilleul              <NA>    CIMETIERE           Arbre
## 9           Cyprès              <NA>    CIMETIERE           Arbre
## 10      Marronnier              <NA> PERIPHERIQUE           Arbre
##    hauteurenm varieteoucultivar         espece
## 1          10              <NA>        cordata
## 2          14              <NA> pseudoplatanus
## 3           0 ''Atropurpureum'' pseudoplatanus
## 4          20              <NA>  hippocastanum
## 5          14              <NA>  hippocastanum
## 6           0              <NA>        pendula
## 7           0              <NA>      americana
## 8          10              <NA>      tomentosa
## 9           0              <NA>      leylandii
## 10          5              <NA>  hippocastanum

Projeter et Afficher les données

La librairie sf surcharge la fonction d’affichage par défaut de R , plot, pour afficher la géométrie des données comme un objet géographique et non comme un nuage de points ou une courbe. Auparavant , nous devons fixer le système de coordonnées de références de l’objet, pour que les données soient correctement projetées à l’affichage.

La fonction st_crs() appliquée sur un objet spatial retourne son CRS s’il est défini , ou permet de fixer avec l’opérateur d’affectation <-. Nous allons utiliser Lambert 93 (EPSG 2154)

st_crs(arbres) 
## Coordinate Reference System:
##   EPSG: 4326 
##   proj4string: "+proj=longlat +datum=WGS84 +no_defs"
arbres <- st_transform(arbres,2154)
arbres <- st_set_crs(arbres, 2154)
st_crs(arbres)  # nouveau CRS de l'objet après transormation
## Coordinate Reference System:
##   EPSG: 2154 
##   proj4string: "+proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

Nous pouvons maintenant afficher les données avec la fonction plot

plot(arbres)

C’est très long ! Celà est dû au fait que la fonction plot affiche autant de cartes que de variables (colonnes) de l’objet spatial.

pour le moment nous n’allons que tracer la géométrie des données : l’attribut geometry de l’objet

plot(arbres$geometry, cex = 0.01)

Tout autre attribut peut être choisi , la fonction de dessin de R se charge de trouver une échelle de couleur en fonction des valeurs que prend la varaible.

plot(arbres["arrondissement"], cex = 0.01, graticule=T)

La sequence de commande suivantes transforme l’attributs arrondissement en facteurs et filtre les données de façon à ne conserver que les arrondissements de Paris intra muros.

arbres$arrondissement <- as.factor(arbres$arrondissement)
levels(arbres$arrondissement) #pour connaître les modalités de la variable
##  [1] "BOIS DE BOULOGNE"  "BOIS DE VINCENNES" "HAUTS-DE-SEINE"   
##  [4] "PARIS 10E ARRDT"   "PARIS 11E ARRDT"   "PARIS 12E ARRDT"  
##  [7] "PARIS 13E ARRDT"   "PARIS 14E ARRDT"   "PARIS 15E ARRDT"  
## [10] "PARIS 16E ARRDT"   "PARIS 17E ARRDT"   "PARIS 18E ARRDT"  
## [13] "PARIS 19E ARRDT"   "PARIS 1ER ARRDT"   "PARIS 20E ARRDT"  
## [16] "PARIS 2E ARRDT"    "PARIS 3E ARRDT"    "PARIS 4E ARRDT"   
## [19] "PARIS 5E ARRDT"    "PARIS 6E ARRDT"    "PARIS 7E ARRDT"   
## [22] "PARIS 8E ARRDT"    "PARIS 9E ARRDT"    "SEINE-SAINT-DENIS"
## [25] "VAL-DE-MARNE"
arbres_intramuros <-  filter(arbres, !(arrondissement %in% c("BOIS DE BOULOGNE",  "BOIS DE VINCENNES", "HAUTS-DE-SEINE", "SEINE-SAINT-DENIS", "VAL-DE-MARNE")))
arbres_intramuros$arrondissement <-  as.character(arbres_intramuros$arrondissement) 
arbres_intramuros$arrondissement <- as.factor(arbres_intramuros$arrondissement)
levels(arbres_intramuros$arrondissement)
##  [1] "PARIS 10E ARRDT" "PARIS 11E ARRDT" "PARIS 12E ARRDT"
##  [4] "PARIS 13E ARRDT" "PARIS 14E ARRDT" "PARIS 15E ARRDT"
##  [7] "PARIS 16E ARRDT" "PARIS 17E ARRDT" "PARIS 18E ARRDT"
## [10] "PARIS 19E ARRDT" "PARIS 1ER ARRDT" "PARIS 20E ARRDT"
## [13] "PARIS 2E ARRDT"  "PARIS 3E ARRDT"  "PARIS 4E ARRDT" 
## [16] "PARIS 5E ARRDT"  "PARIS 6E ARRDT"  "PARIS 7E ARRDT" 
## [19] "PARIS 8E ARRDT"  "PARIS 9E ARRDT"

On peut maintenant tracer les positions des arbres de Paris intra-muros, après avoir défini le CRS de cette couche par la fonction st_crs()

st_crs(arbres_intramuros)
## Coordinate Reference System:
##   EPSG: 2154 
##   proj4string: "+proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
plot(arbres_intramuros["arrondissement"], cex=0.1)

Données vectorielles des quartiers de Paris.

Nous allons maintenant charger les contours des quartiers de Paris, disponibles sur (https://opendata.paris.fr/explore/dataset/quartier_paris/information/). Chaque arrondissement contient 4 quartiers, on a donc 80 unités spatiales à considérer.

quartiers <-(read_sf("quartier_paris.shp"))
quartiers <- st_transform(quartiers, 2154)
quartiers <-  st_set_crs(quartiers,2154)

plot(st_geometry(quartiers))
plot(arbres_intramuros["arrondissement"], add=TRUE, alpha=0.5, cex=0.1  )

Cartographie simple

Pour tracer un objet spatial , il suffit d’appliquer la fonction plot() sur cet objet. Cette fonction peut également colorer la géométrie de l’objet en fonction d’une variable

Les 6 genres d’arbres les plus représentés

Nous utilisons les fonctions countdu package dplyr, pour compter les nombre d’arbre par genre. (c’est l’quivalent d’un group_by, suivi d’une aggregation count) Nous utilisons la fonction top_n du package dplyr pour selectionner les 6 genres les plus représentés. Nous utilisons la fonction filter et l’operateur %in% pour ne conserver que les arbres de ces six genres.

count_by_genre <-  count(arbres_intramuros,genre, libellefrancais, sort = T,name = "nb_indiv")
head(count_by_genre, 10)
## Simple feature collection with 10 features and 3 fields
## geometry type:  MULTIPOINT
## dimension:      XY
## bbox:           xmin: 644413.8 ymin: 6857489 xmax: 657170.6 ymax: 6867050
## epsg (SRID):    2154
## proj4string:    +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## # A tibble: 10 x 4
##    genre  libellefrancais  nb_indiv                                geometry
##    <chr>  <chr>               <int>                        <MULTIPOINT [m]>
##  1 Plata… Platane             37856 (645035 6861100, 645068.5 6861382, 645…
##  2 Aescu… Marronnier          20051 (644413.8 6861116, 644420.5 6861114, 6…
##  3 Tilia  Tilleul             16778 (645059.1 6860934, 645070.1 6861108, 6…
##  4 Acer   Erable              12369 (645032.2 6860884, 645037.9 6860896, 6…
##  5 Sopho… Sophora             10708 (645069.1 6861372, 645088.8 6860152, 6…
##  6 Pinus  Pin                  3751 (645037.9 6860889, 645057.3 6861095, 6…
##  7 Celtis Micocoulier          3628 (645102.9 6860998, 645167.5 6861102, 6…
##  8 Fraxi… Frêne                3580 (645068.1 6860577, 645074.4 6861408, 6…
##  9 Prunus Cerisier à fleu…     3359 (645037.9 6860884, 645056.4 6860576, 6…
## 10 Pyrus  Poirier à fleurs     3230 (645093.6 6861056, 645103.2 6861005, 6…
six_genres <-  top_n(count_by_genre,6,nb_indiv)
head(six_genres)
## Simple feature collection with 6 features and 3 fields
## geometry type:  MULTIPOINT
## dimension:      XY
## bbox:           xmin: 644413.8 ymin: 6857489 xmax: 657126.1 ymax: 6867050
## epsg (SRID):    2154
## proj4string:    +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## # A tibble: 6 x 4
##   genre  libellefrancais nb_indiv                                  geometry
##   <chr>  <chr>              <int>                          <MULTIPOINT [m]>
## 1 Plata… Platane            37856 (645035 6861100, 645068.5 6861382, 64509…
## 2 Aescu… Marronnier         20051 (644413.8 6861116, 644420.5 6861114, 644…
## 3 Tilia  Tilleul            16778 (645059.1 6860934, 645070.1 6861108, 645…
## 4 Acer   Erable             12369 (645032.2 6860884, 645037.9 6860896, 645…
## 5 Sopho… Sophora            10708 (645069.1 6861372, 645088.8 6860152, 645…
## 6 Pinus  Pin                 3751 (645037.9 6860889, 645057.3 6861095, 645…

L’immense avantage du package sf est que que les géométries ont été conservées lors des manipulation de regrouppement , filtrage etc. Il n’est donc pas nécessaire de refaire les opérations sur les géométries , qui ont été transmises aux résultats lors des opérations count et top_n.

plot(six_genres["libellefrancais"], 
     key.pos=1, cex=0.2, 
     key.length = 0.7,
     main="Les six genres d'arbres les plus représentés à Paris")

Que peut-on dire de l’implantation de ces 6 genres ?

La domanialité

Afficher cette variable est immédiat. La légende est perfectible pour le moment; pour faire une vraie carte , voir à la fin du support l’utilisation du package cartography.

plot(arbres_intramuros["domanialite"], 
     key.pos=4, cex=0.2, 
     key.length = 1,
     key.width = lcm(4.5),
     main="Domanialité des arbres de Paris")

Calcul du nombre d’arbres par quartier

Aggregation spatiale

On peut faire la jointure spatiale avec le dataframe arbres_intramuros, pour inclure les informations de quartiers aux arbres à l’aide de la relation spatiale st_within

arbres_in_quartiers <- st_join(arbres_intramuros,quartiers, join=st_within)
nrow(arbres_in_quartiers) 
## [1] 163367

On obtient alors un nouveau dataframe spatial, contenant les 163367 arbres situés dans les quartiers intramuros, et augmentés des variables issues du dataframe quartier.

On peut alors compter le nombre d’arbres par quartiers avec la fonction table qui compte (entre autres) le nombre d’individus (ligne) par valeurs distinctes d’une de ses variables. On peut étendre le nombre de variables pour obtenir des tables de contingence.

nb_arbres_by_quartier <- table(arbres_in_quartiers$c_qu)

Pour ajouter cette information au dataframe quartiers, l’une des possibilité est de trier le dataframe par son code de quartier c_qu puis d’ajouter la colonne nb_arbres_by_quartiers obtenue précédemment.

on peut par exemple utiliser la fonction arrange du package dplyr

quartiers <-  arrange(quartiers,c_qu)
quartiers$nb_arbres <-  nb_arbres_by_quartier

Avant de cartographier cette valeur, regardons sa distribution à l’aide d’un histogramme pour déterminer si un traitemment particulier doit être envisagé (en cas de distribution particulièrement biscornue)

hist(quartiers$nb_arbres,breaks = 10)

Carte choroplète du nombre d’arbres (mauvaise pratique)

plot(quartiers["nb_arbres"], main=NULL, key.pos = 4)
title("Nombre d'arbres par quartier administrif de Paris",sub = "Représentation sémiologiquement discutable" )

“Raster” du nombre d’arbres

Pour réaliser une sorte de raster, nous allons créer une grille qui couvre la zone d’étude, et projeter les points du semis (le dataframe arbres_intramuros) dans les cellules créées (c’est également une façon d’approcher une carte de densité 2D de façon discrète)

Afin de ne pas faire une grille qui couvre la boîte englobante de la zone d’étude, mais que la grille ne couvre que les géométries de chaque quartier, on réalise une intersection entre la grille raster et l’enveloppe des quartiers.

N.B. ce n’est pas la façon canonique de réaliser un raster, pour cela il faut utiliser la library raster et réaliser le raster à partir l’aide de la fonction

rast1 <- st_make_grid(quartiers, square = TRUE, n=40, what="polygons") ## raster de 1600 cellules
plot(quartiers$geometry)
plot(rast1,  add=T)

#restriction du raster à la geométrie des quartiers par intersection avec son enveloppe
envelop_qu <-  st_union(quartiers)
plot(envelop_qu)

rast2 <-  st_intersection(envelop_qu, rast1) ## raster de 1065 cellules

# on retransforme l'objet rast2 qui est une collection (sfc) en objet sf
rast2 <- st_sf(rast2)
plot(rast2)

#on récupère les points contenus dans les cellules
predicat_intersect <- st_contains(rast2,arbres_intramuros)

#on compte le nombre de points par cellules avec la fonction length
rast2$nb_arbres <- unlist(lapply(predicat_intersect, length)) 
hist(rast2$nb_arbres)

plot(rast2)

A quoi correspondent les cellules jaunes ?

On peut réitérer le processus en changeant la taille de la grille ou le type de grille (hexagonale)

rast1 <- st_make_grid(quartiers, square = TRUE, n=80)
rast2 <-  st_intersection(envelop_qu, rast1) 
rast2 <- st_sf(rast2)
predicat_intersect <- st_contains(rast2,arbres_intramuros)
rast2$nb_arbres <- unlist(lapply(predicat_intersect, length)) 

plot(quartiers$geometry, border="white", bgc="#222222"  )
plot(rast2, border=NA, key.pos=4, add=T)
plot(quartiers$geometry, border="white", bgc="#222222", add=T, alpha=0.8,lwd=0.3,key.pos=4)

rasthex <- st_make_grid(quartiers, square = FALSE, n=60, what="polygons")
rasthex <- st_sf(rasthex)
#les hexagones s'arrètent à l'evloppe, pas besoin de la créer
predicat_intersect <- st_contains(rasthex,arbres_intramuros)
rasthex$nb_arbres <- unlist(lapply(predicat_intersect, length)) 

plot(quartiers$geometry,bgc="#222222")
plot(rasthex,  key.pos=4, add=T, lwd=0.2)
plot(quartiers$geometry, border="white", bgc="#222222", add=T, alpha=0.8,lwd=0.3,key.pos=4)

On peut noter que le changement de découpage de l’espace dû à la taille ou la forme des cellules change l’aspect de la carte.

Calcul et cartographie de la densité d’arbres par quartiers

Pour calculer la densité d’arbres par quartier, nous devons compter le nombre d’arbres par quartier, et diviser ce nombre par la surface du qartier.

Nous avons déjà ajouté un attribut nb_arbres à l’objet quartiers. L’objet quartier contient aussi un attribut surface, issu des données brutes. On doit d’abord vérifier que cette valeur est correcte (la projection des données est équivalente) en calculant l’aires des quartiers et en les comparant à l’attribut pré-existant

ecarts <-  as.numeric(st_area(quartiers)) - quartiers$surface
ecarts
##  [1] 1.877930e-05 8.108851e-06 6.117160e-06 5.210517e-06 4.029163e-06
##  [6] 5.435315e-06 6.035494e-06 6.569724e-06 7.249997e-06 5.676586e-06
## [11] 8.165312e-06 4.341156e-06 6.586371e-06 9.512180e-06 1.055695e-05
## [16] 8.231029e-06 1.244282e-05 1.792773e-05 1.491560e-05 9.088253e-06
## [21] 6.806862e-06 1.615961e-05 1.833285e-05 5.744630e-06 1.758547e-05
## [26] 2.283882e-05 1.835986e-05 2.939277e-05 2.566329e-05 1.731119e-05
## [31] 1.683447e-05 2.600718e-05 1.549569e-05 1.163094e-05 8.670206e-06
## [36] 1.081359e-05 2.129877e-05 9.055191e-06 1.324248e-05 1.921703e-05
## [41] 1.531083e-05 1.852063e-05 2.501090e-05 2.072949e-05 1.291232e-04
## [46] 1.569642e-04 4.206994e-05 2.605096e-05 2.508820e-05 6.487360e-05
## [51] 4.781876e-05 1.511956e-05 2.415944e-05 2.820021e-05 2.937065e-05
## [56] 3.764336e-05 6.029010e-05 3.513182e-05 3.172900e-05 5.560508e-05
## [61] 1.362823e-04 1.181327e-04 6.534485e-05 3.181025e-05 3.170897e-05
## [66] 2.902956e-05 3.089476e-05 2.964167e-05 4.153443e-05 3.538677e-05
## [71] 2.405327e-05 2.895016e-05 2.844469e-05 5.151751e-05 3.934582e-05
## [76] 2.758135e-05 1.779071e-05 3.266637e-05 3.387243e-05 4.496775e-05

Les écarts sont très faibles (de l’ordre de \(10^{-5}m^2\)) , suffisamment pour que l’attribut surface soit directement utilisé pour le calcul de la densité.

On crée une variable dens dans l’objet quartier dont la valeur est le ratio entre la variable surface et la variable nb_arbres. On aura vérfié auparavant qu’il n’y a pas de valeurs manquantes dans l’objet , avec la fonction anyNA()

quartiers$dens <- quartiers$nb_arbres / quartiers$surface
plot(quartiers["dens"], main="Densité d'arbres par quartier")

Connaissant l’implantation des arbres et la formes des quartiers , que dire des densité des quartiers sud-ouest et sud-est ?

Quel est le quartier le plus dense ?

quartiers %>% top_n(1,dens)
## Simple feature collection with 1 feature and 10 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 654747.9 ymin: 6862008 xmax: 656445.5 ymax: 6863585
## epsg (SRID):    2154
## proj4string:    +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## # A tibble: 1 x 11
##   n_sq_qu  c_qu c_quinsee l_qu   c_ar n_sq_ar perimetre surface
## *   <dbl> <dbl>     <dbl> <chr> <dbl>   <dbl>     <dbl>   <dbl>
## 1  7.50e8    79   7512003 Père…    20  7.50e8     5014.  1.60e6
## # … with 3 more variables: geometry <POLYGON [m]>, nb_arbres <int>,
## #   dens <dbl>

Cartographie par symboles proportionnels

Pour utiliser des symboles proportionnels , nous devons utiliser un package spécifique : cartography. Dans ce package , l’ajout de couche au dessin courant est activé par défaut (add=TRUE)

library(cartography)
plot(quartiers$geometry, bgc="#888888", border="white", lwd=0.3)
propSymbolsLayer(quartiers,var = "nb_arbres", inches = 0.15, legend.pos = NA)
# couches de disposition des légendes et du texte
layoutLayer(title = "Nombre d'arbres à Paris par quartier",  
            frame =F, north = TRUE, author="M2 IGAST2019-2020",
            sources = "Opendata.paris.fr 2019", 
            scale = 50)
legendCirclesSymbols(pos = "topleft", title.txt = "nombre d'arbres", title.cex = 0.8, cex = 1,
  border = "black", lwd = 1, values.cex = 0.6, var=c(50,1000,3500,7000), inches=0.15,
  col = "#E84923", frame = FALSE, values.rnd = 0, style = "c")

Le package cartography est un excellent package, que je recommande vivement.